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THE  MULTIVARIATE  GRAM-CHARLIER  SERIES  APPLIED  TO 
RANDOM  SIGNAL  DETECTION 

1.  INTRODUCTION 


The  application  of  probability  theory  and  statistical  modeling  to  the  analysis  of  physical 
systems  frequently  leads  to  mathematical  problems  that  are  not  amenable  to  closed-form  solution. 
Complicated  nonlinearities,  unwieldy  computations,  and  other  manipulative  difficulties  often 
necessitate  the  use  of  numerical  techniques  and/or  approximations.  In  the  latter  case,  several 
alternatives  are  possible  One  common  approach  is  to  represent  the  probability  density  function 
(PDF)  as  an  infinite  series  of  orthogonal  polynomials,  with  the  expansion  coefficients  being 
specified  in  terms  of  statistical  moments. The  key  to  successful  application  of  this  procedure 
lies  in  judiciously  selecting  the  orthogonal  polynomials.  For  practical  applications,  series 
approximations  prove  useful  only  when  a  small  number  of  terms  suffice  to  provide  the  required 
accuracy.  In  this  context,  the  Gram-Charlier  series^  is  known  to  render  good  service,  especially  in 
situations  where  the  original  PDF  exhibits  Gaussian  characteristics;  viz.,  unimodality,  continuity, 
and  bounded  variation.  The  theoretical  and  computational  aspects  of  approximating  a  univariate 
PDF  via  the  Gram-Charlier  series  have  been  extensively  documented  in  the  literature.'’^ 
Unfortunately,  corresponding  results  for  a  multivariate  PDF  are  conspicuously  meager.  This  is 
due  to  the  added  complexity  that  typically  characterizes  «-dimensional  mathematical  analysis,  as 
well  as  a  general  lack  of  familiarity  with  the  properties  of  orthogonal  polynomials  in  several 
variables.  Such  considerations  need  not  be  of  great  concern,  however,  since  the  Gram-Charlier 
series  expansion  technique  can  actually  be  extended  to  higher  dimensions  without  explicit  use  of 
complicated  analytical  procedures  or  functions;  viz.,  tensor  analysis^  and  multivariate  orthogonal 
polynomials.’  A  straightforward  extension  is  described  in  the  first  part  of  this  report. 

Section  2. 1  focuses  on  developing  a  suitable  infinite  series  representation  for  the  n- 
dimensional  Dirac  delta  function  The  Gram-Charlier  expansion  of  an  arbitrary  multivariate  PDF 
is  then  obtained  as  demonstrated  in  section  2.2.  There,  analytical  expressions  for  the  unknown 
expansion  coefficients  are  also  developed  and  explicitly  evaluated  in  terms  of  the  central  moments 
of  the  distribution.  In  section  2.3,  a  linear  integral  equation  of  the  second  kind  is  derived  for 
subsequent  use  in  solving  an  important  binary  detection  problem  that  arises  in  sonar  array 
processing  applications.  Finally,  section  2.4  provides  a  comprehensive  discussion  of  the 
mathematical  results 

The  second  part  of  this  report  examines  the  problem  of  detecting  random  multivariate  signals 
embedded  in  Gaussian  noise.  A  canonical  representation  for  the  likelihood  ratio  is  derived  in  the 
form  of  an  infinite  series  whose  terms  depend  on  the  received  measurement  vector,  central 
moments  of  the  random  signal,  and  the  background  noise  statistics.  The  representation  explicitly 
prescribes  the  optimum  processing  required  to  detect  a  broad  class  of  multivariate  random  signals. 
Unfortunately,  few  representations  of  the  type  developed  here  are  available  in  the  literature. 

Those  that  do  appear  usually  avoid  detailed  characterization  of  non-Gaussian  signals.  Low  signal- 
to-noise  ratio  conditions  are  typically  assumed  to  exist,  allowing  the  use  of  simplifying 
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approximations  that  preclude  explicit  specification  of  the  signal  PDF;  only  its  mean  and 
covariance  are  required,  as  in  the  case  of  locally  optimal  Bayes  detectors.*’^  A  notable  exception 
to  the  aforementioned  approach  is  the  work  of  Kelly  and  Tague.*®  There,  a  canonical 
representation  is  developed  for  non-Gaussian  univariate  signals  without  any  restrictions  on  input 
signal-to-noise  ratio.  The  only  assumptions  made  are  that  the  background  noise  is  Gaussian 
distributed  and  the  signal  PDF  can  be  represented  by  a  Gram-Charlier  series  expansion.  The  first 
assumption  is  reasonable  for  most  sonar  signal  processing  applications,  since  ambient  noise  in  the 
ocean  environment  usually  exhibits  Gaussian  statistical  characteristics  The  second  assumption 
imposes  no  significant  restrictions  either,  since  it  is  satisfied  by  a  wide  variety  of  acoustic  signals 
encountered  in  practice.  In  fact,  the  generality  of  Kelly  and  Tague's  work  renders  it  an  important 
contribution  to  the  signal  processing  literature. 

Despite  its  obvious  attributes,  the  aforementioned  representation  was  developed  for 
univariate  signals  only;  thus  it  cannot  be  used  effectively  for  array  processing  applications 
involving  vector-valued  measurements.  A  desire  to  overcome  this  limitation  provides  impetus  for 
the  present  study.  Following  the  approach  originally  employed  by  Kelly  and  Tague,'®  pertinent 
mathematical  formulae  are  developed,  extended,  and  modified  to  accommodate  multivariate 
signals.  The  resulting  n-dimensional  representation  retains  all  the  functional  characteristics, 
attributes,  and  generality  of  its  one-dimensional  counterpart,  and  as  expected,  becomes  identical 
to  it  for  the  case  n  =  1 .  An  attractive  feature  of  all  formulae  is  that  statistical  means  and 
covariances  associated  with  signal  and  noise  always  appear  in  their  compact  vector  and  matrix 
formats,  respectively.  The  need  to  explicitly  depict  individual  vector  components  or  matrix 
elements  is  avoided  via  utilization  of  Kronecker  products.'^  Indeed,  this  artifice  greatly  simplifies 
expressions  that  would  otherwise  be  cumbersome  and  not  amenable  to  analysis  and/or 
interpretation. 

Section  3.1  provides  a  mathematical  formulation  of  the  binary  detection  problem  for  random 
multivariate  signals  embedded  in  Gaussian  background  noise.  The  problem  is  examined  in  the 
context  of  sonar  array  signal  processing,  and  the  solution  is  specified  as  a  statistical  hypothesis 
test;  viz.,  the  "likelihood  ratio  test."^’'°’'^  In  section  3.2,  a  canonical  representation  for  the 
likelihood  ratio  is  developed.  This  representation  takes  the  form  of  an  infinite  series  whose  terms 
depend  on  the  received  measurement  vector  together  with  the  signal  and  noise  statistics.  Some 
interesting  possibilities  for  approximating  the  likelihood  ratio  are  discussed  in  section  3.3.  The 
various  alternatives  are  revealed  via  analysis  of  the  series  expansion.  In  addition,  insight  is  gained 
regarding  the  optimal  signal  processing  required  to  detect  non-Gaussian  signals  embedded  in 
Gaussian  noise. 

It  is  not  necessary  to  go  through  the  mathematical  derivations  presented  in  the  first  part  of 
this  report  completely  before  proceeding  to  the  signal  detection  material  discussed  in  the  second 
part.  The  reader  interested  primarily  in  results  and  applications  may  omit  the  first  part  entirely,  or 
content  himself  with  making  only  slight  acquaintance  with  the  multivariate  Gram-Charlier  series 
derivation.  Pertinent  results  used  in  the  second  part  are  appropriately  referenced  to  the  first  part 
as  they  occur. 
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2.  THE  MULTIVARIATE  GRAM-CHARLIER  SERIES 


2.1  REPRESENTATION  OF  THE  /i-DIMENSIONAL  DIRAC  DELTA  FUNCTION 
It  is  well  known**'  that  the  functions 

P,  (z)  =  (2*  ^  H,  (z)  k  =  0,1,2 

where  H,^  (z)  are  Hermite  polynomials  defined  by  the  formulae*** 


A:  =  0,1,2,..., 


form  an  orthonormal  system  on  the  interval  -00  <  z  <  00 ,  Accordingly,  these  functions  must 
satisfy  Bessel's  equality,*^  or  the  equivalent  completeness  relation 


2-*r  ^ 


where  S[z  -  z)  represents  the  one-dimensional  Dirac  delta  function  * 


If  one  makes  the  substitutions 


u  =  u.  /  V2 


for  u  =  (?,z), 


in  equation  (3),  and  uses  the  chain  rule  formula  for  derivatives 


^  ^  =  J2  ^ 


forM  =  (z,z), 


along  with  the  relation* 


'  Throughout  this  report,  the  notation  employed  in  equation  (4)  will  be  used.  The  desired  substitutions  are 
z  =  z"  /  V2  and  z  =  z,  /  42  .  To  obtain  these  fomulae,  simply  replace  the  letter  u  with  z  and  z,  respectively. 


it  is  readily  shown  that 


S(z,-z,)=^  -y-L 
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For  notational  convenience,  one  can  introduce  the  scalar  function 


for  u  -  {z  ,z), 


and  the  differential  operator 


These  quantities  allow  equation  (7)  to  be  recast  in  a  more  compact  form; 

- ^. )  =  S  A ^ ^ 

k=o  ^  ’ 

Although  equations  (3)  and  (10)  describe  identical  completeness  relations  in  one  dimension, 
the  latter  expression  is  more  readily  extended  to  «-dimensions.  To  this  end,  let 


for  M  =  {z,z). 


where  the  superscript 'T'  denotes  matrix  transpose,  and  observe  that 


-]/2(u'  u) 


for  u  =  {z,z), 


n^(z;-z,)  =  (5(z-z). 


Using  the  preceding  formulae  in  conjunction  with  equation  (10)  leads  to  the  result 


Wz)o(z)  j(z-z)=  n  zw 

,=1  \ki=0  •  J 
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The  product  of  n  infinite  series  appearing  in  square  brackets  on  the  right-hand  side  of 
equation  (14)  can  be  expressed  as  a  single  infinite  series  by  repeated  application  of  the  double 
summation  formula'® 


ZA  is. 

k~0  y  \m=0  '  m=0/c=0 


To  see  this,  let 


A=n  s 


^  ^  k  I 

i=pyk,=o 


p  =  \,2,...,n. 


Equation  (14)  then  takes  the  form 

-  z)  =  /,<D(z  )0(z).  ( 1 7) 

It  now  follows  from  equation  (15),  equation  set  (16),  and  the  Binomial  Expansion  Theorem'®  that 


(a)"Y^(4)' 


k  ^(k  —  ^  V  ^ 
ki=Q  -V^2  J 


kj  =  0  ^2'  k^=0  '\^2  ^\)  J 


I  h- 


^,=0  K'- 


Repeating  this  procedure  yields  a  similar  expression 
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f  oo  1  M  ^  ( T 

^2  J\k^=0  ^3-  ) 


Vkh  k,!{k,-k,)i  J- 


/cj=0  ^3-  Jt2=0  ^2-V^3  ^2/  J 


X  — (i,  +  i,+4)‘  /. 


^  ^  f 

k.^O  ^3  • 


and,  finally,  by  induction  the  following  result  is  obtained: 


The  desired  single  series  representation  for  /i  is  now  deduced  by  noting  that 


yL=y  ^.■■.=[-^1  \  — 

^  ^  \(9FJ  \<^Zy 


where 


A-  A  A.  A 

du  du^’  du.^'  '  du„ 


for  M  =  {z,z), 


is  the  gradient  operator  associated  with  the  vector  u.  Combining  the  results  embodied  in 
equations  (17),  (20),  and  (21)  yields  the  ^-dimensional  analogue  of  equation  (10): 


[£]  [fj 


or,  equivalently. 
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f  SY  f 

5{l-z)  -A  1 
/(D(z)<D(z)  ”§A' 


To  facilitate  the  Gram-Charlier  expansion  of  an  arbitrary  multivariate  PDF  later  on,  it  is 
convenient  at  this  point  to  change  the  variables  from  z  and  z  to  x  and  x  via  the  relations 

z  =  r'^'(x-x), 


z  =  r'^^(x-x),  (26) 

where  F  is  an  arbitrary  («  x  n)  symmetric,  positive  definite  matrix,  and  x  is  an  arbitrary  («  x  1) 
vector  (these  quantities  will  be  specified  later).  Substituting  equations  (25)  and  (26)  into  equation 
(24),  and  using  the  following  relationships 


S{r^'^{x-x})=^\S{x-x), 


d  ^pi/2  d 
dz  dx 


r  =  r'/T'/\ 


leads  to  the  expression 


^(x  -  x) 


1 IIJ  <3 


A^(x;x.r)A^(x;x,r) 


A^(x;x,r)A^(x;x,r) 


where 
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iv(w,x,r)  = 


for  u  =  (x,x), 


is  the  w-dimensional  Gaussian  PDF, 

Although  the  series  representation  for  ^(x  -  x)  given  by  equation  (32)  is  pleasingly  compact, 

it  does  not  yet  exhibit  the  "separation  attribute"  necessary  for  utilitarian  applications.  Specifically, 
individual  terms  on  the  right-hand  side  of  equation  (32)  are  not  written  in  a  form  that  explicitly 
uncouples  the  x  dependence  from  the  x  dependence.  This  can  be  accomplished,  however,  by 
using  the  Kronecker  product  formula  for  differential  operators  derived  in  appendix  A  (see 
equation  (A- 17));  viz.. 


— 1  r(  — 

<7xJ  v<7x 


11:  =  0,1,2,... 


where  the  notation 


[yf®...]*  =J<E>A0...<8>A 


yt  =  0,1,2,..., 


has  been  employed.  Substituting  equation  set  (34)  into  equation  (32)  yields  the  desired  result 

iV(x.-x.r)A(x.-x.r) 


Here, 


N{u;x,T) 


^  =  0,1,2,..., 


for  M  =  (x ,  x). 


is  a  («^  X  1)  vector  valued  function  of  its  argument  u.  An  alternative  representation  for 
^^[u;x,T),  which  is  simpler  to  evaluate  explicitly,  may  also  be  derived  by  combining  equations 

(12),  and  (25),  through  (31)  with  equations  (33)  and  equation  set  (37).  Performing  the  necessary 
algebraic  manipulations  yields 


I,  (w,-  X,  r)  =  ^^  (r-’"'  {u  -  x))  k  =  0,1,2, . . . 


for  u  =  (x,x), 
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where 


]  f 

y;  (u)  =  -—  —  o...  0(m) 

0(w)  Vdu) 


k  =  0,1,2, 


2.2  EXPANSION  OF  AN  ARBITRARY  MULTIVARIATE  PDF 

Let  X  be  a  random  («  x  1)  vector  defined  over  all  A7-dimensional  space  R„  whose  PDF  is  given 
by  p(^.  The  mean  value  of  x  and  its  covariance  matrix  then  satisfy  the  relations 

X  -  f  ^(x)afx,  (40) 

r  =  (x  -  x)(x  -  x)^  p(x)^2E-  (4 1 ) 

Substituting  these  parameters  into  equation  (36),  multiplying  both  sides  of  that  expression  by 
p  (x ) ,  and  then  integrating  the  result  over  Rn  (here  x  is  used  as  the  variable  of  integration)  yields 
the  formula 

p(x)  =  ^^{x;x,T),  (42) 

A:^0  ^  • 

where 

Equation  (42)  describes  the  Gram-Charlier  series  expansion  of  an  arbitrary  multivariate  PDF. 
The  expansion  coefficients  are  specified  in  terms  of  the  central  moments  of  the  PDF  via 

equation  (43).  To  compute  P^  explicitly,  it  is  first  necessary  to  determine  the  vector-valued 

polynominals*  i^^(x,  x,r)  in  closed  form.  To  this  end,  recall  that  ^^(x,  x,r)  is  simply  related  to 

V'tlr  '  ^  (x-x})  as  is  indicated  by  equation  set  (38).  Consequently,  knowledge  of  the  latter 
provides  sufficient  information  to  evaluate  the  former.  Since  ^ satisfies  equation  set  (39),  it 
immediately  follows  that 


The  components  of  ^  (x,  x,  T )  are  related  to  the  Hermite  polynomials  in  several  variables  described  by 


Erdelyi. 
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1  d 

- ( 


...  o(m) 


1  d 

- 1 

3)(w)  du 


)[<D(wV^(m)] 


<I)(m)  du  — *  <?w  ~'‘ 


=  -w0^^(m)  +— (8)^^(m). 

Using  the  functional  definition  embodied  in  equation  set  (39),  the  recursion  formula  specified  by 
equation  (44),  and  the  differentation  rules  given  in  appendix  A,  it  is  possible  to  determine 

polynomial  representations  for  all  members  of  the  set  A:  =  0,1,2,...|.  Note,  however,  that 

If/ ^  (m)  is  a  («*"  X  1)  vector  valued  function  of  its  argument  u.  Unfortunately,  since  the 

dimensionality  of  grows  exponentially  with  increasing  index  k,  there  is  a  concomitant 

increase  in  computational  complexity  associated  with  the  explicit  evaluation  of  each  successive 
vector.  To  forestall  this  encumbrance,  an  alternative  recursion  formula  may  be  employed.  In 
particular,  observe  that  (see  appendix  A) 


0(m)  = 


du\du 


0(m)  [v  du. 


'•  •  •  ^  (w)  ■ 


0(m)  du 


—  —  ®...  0(m) 

duVduJ 


where 


duydu 


/,y  =  1.2,...,«, 
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is  a  («  X  ri)  matrix  operator  whose  {i,j)  element  is - .  This  matrix  is  sometimes  referred  to  in 


the  literature  ’  as  the  "Hessian"  operator.  If 


1  ^  f  ^v) 

0(m)  J 


^2k+l(—)  ~ 


1  ^ 

- I 

<[)(«) 


AfA 

du  V  du 


5...  a)(M) . 


it  then  follows  from  equation  set  (39)  and  equations  (45),  (46),  (48),  and  (49)  that 


yt  =  0,1,2,... 


A  recursive  relation  for  derives  by  substituting  equation  (48)  into  equation  (49). 

The  result  is 


Next,  the  matrix  transpose  of  equation  (49)  is  taken  and  the  Hessian  operator  symmetry  property 


^  f  d 


du\du 


A\A 

du  V  du 


is  used  to  obtain  the  expression 


(53) 


(55) 


Finally,  letting  A:  =  0  in  equation  (48)  leads  to  the  initial  condition 

M„(w)=l.  (56) 

Equations  (51),  (55),  and  (56)  specify  a  recursive  procedure  for  computing  the  set  of 
matrices  (m):  k  =  0,1,2,. .  .  An  examination  of  equations  (48)  and  (49)  reveals  that  the 

dimension  of  M2*(m)  and  A/2t+i(w)  are  of  dimension  («*  x  «*)  and  x  w*).  Furthermore,  since 
the  vec  (•)  operator  simply  rearranges  matrix  elements  in  a  prescribed  fashion  (see  appendix  A), 
it  follows  from  equation  set  (50)  that  can  be  easily  evaluated  once  Mk(u)  is  determined 

explicitly.  At  the  2k'^  stage,  actual  computations  involve  (w*"  x  «*)  dimensional  matrices  rather 
than  (a?^  X  1)  dimensional  vectors.  This  reduction  in  the  exponential  growth  rate  by  a  factor  of  2 
significantly  reduces  the  complexity  of  required  mathematical  operations.  For  convenience,  the 
pertinent  equations  and  first  few  polynomial  representations  are  summarized  as  follows: 
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(57a) 

(m)  =  -  w  ®  («)  +  ;^  ®  ^2*  («). 

(57b) 

H..2  (w)  =  -H  ®^2lM  +  •£  ®'”2L(")- 

(57c) 

y/ J(^  =  vec^M^{i^  ^  =  0,1,2,..., 

(57d) 

M,(«)  =  1, 

(58a) 

M,  (m)  =  -M, 

(58b) 

M2{u)  =  ut/  -  I, 

(58c) 

M^(u)  =  C[/0w]- 

(58d) 

K(“)  =  ®A^2(h)]C-C. 

(58e) 

Here, 

I  =Z  £.(")/(«) 

/-I 

(59) 

is  the  (n  X  «)  identity  matrix.  In  addition. 

C  =  I®I  +  U{n,n), 

(60) 

where 

=  S  Z  {^. W) ®  {£>('')£•  W) 

j=l  ;=1  ^ 

(61) 

is  the  {n^  x  «^)  permutation  matrix,*^  and 
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are  (w  x  1)  unit  vectors. 

Returning  to  the  problem  at  hand,  a  closed-form  expression  for  the  vector  valued 
polynomial  follows  immediately  by  combining  equation  set  (38)  with  equation  set 


(50).  In  particular, 

|,(^.|.r)  =  vec{M,(r'"'{x-x})}  A -0,1,2 .  (63) 

Substituting  this  expression  into  equation  (43)  then  yields 

=  vecjj^  -x})/7(xyx|  A:  =  0,1,2,....  (64) 

Use  of  equations  (40),  (41),  and  equation  set  (58)  allows  explicit  evaluation  of  as  shown 
below: 

_  p  =\,  (65a) 

LL.O 

=  0  =  («  X  \)null  vector,  (65b) 

P^  =  000  =  X  ^null  vector,  (65c) 

p^  =  -[r-'^‘0. . .]'£  [(x  - x)0. . .fp{x)dx.  (65d) 


p^  =  [r  "' 0. . .]'  £  [(x  -  x)0. . .]V(x)4/x  -  vec{C(r  0  r)}  -  vec(r) 0  vcc(r)  .  (65e) 

Since  each  term  in  the  multivariate  Gram-Charlier  series  expansion  of /7(x)  takes  the  form  of 
a  scalar  product  of  two  vectors,  the  computation  process  can  be  simplified  further  by  employing 
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various  mathematical  artifices.  Details  of  the  procedure,  along  with  explicit  evaluation  of  the  first 
two  non-trivial  terms,  are  presented  in  appendix  B. 


2.3  AN  INTEGRAL  EQUATION  SATISFIED  BY  THE  MULTIVARIATE 
GRAM-CHARLIER  EXPANSION  FUNCTIONS 


In  section  2. 1,  the  Gram-Charlier  expansion  functions  were  defined  by  the  differential 
relationships  expressed  in  equation  set  (37),  or  equivalently,  in  equation  set  (63).  An  alternative 
relationship  in  integral  form  may  also  be  developed  by  using  the  well-known  characteristic 
function  formula  for  a  n-dimensional  normalized  Gaussian  PDF,^  viz.. 


*(»)= 


e  ~  '  ^ 


(66) 


where  <!)(«)  is  specified  by  equation  (12).  Differentiating  this  expression  successively  with 
respect  to  u,  and  noting  that 


3  T  ^ 

—  ^  -  yW-" 

du 


(67) 


leads  to  the  expression 


r 


\duj 


(2;r)"-’«”“  •' 


(68) 


An  integral  formula  for  y/  ^  (m)  is  now  obtained  by  dividing  both  sides  of  equation  (68)  with 
<I>(m),  and  then  using  equations  (12)  and  (43).  The  result  is 


[Itt) 


(69) 


Combining  equation  set  (38)  with  equation  (69)  also  yields 


(70) 


To  derive  the  integral  equation  for  multiply  both  sides  of  equation  (69)  by 

n(u;u,[1  +  Q]"' j  and  integrate  the  resulting  expression  over  R„to  obtain 
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(71) 


(2a-)"''  (2a)”'' 


(C 


[w0...] 


|i  ->/2|(!;-u)''[/  +  Q]((i-£)-(u+>w)'^(!(+>!S')| 


dwdu. 


Here,  it  is  assumed  that  Q  is  a  positive  definite  symmetric  matrix.  This  assumption  insures  that 
the  integrals  appearing  on  the  right-hand  side  of  equation  (71)  are  both  absolutely  convergent, 
thus  permitting  the  order  of  integration  to  be  interchanged.  Performing  this  operation  allows 
equation  (70)  to  be  rewritten  in  the  form 


A^^w;g,[/ -i-Q]  ^^y/^{u)du 


ill) 


du  uw 


where 


u  =  m  +  Q'(m-i-  yw). 


(73) 


Since  the  integral  enclosed  by  the  square  brackets  on  the  right-hand  side  of  equation  (72)  is  equal 
to  unity,  it  follows  that 


l/2{u  +  y  w)^|/  +  Q"*  j(u+y  w) 


dw,  (74) 


Under  the  assumption  that  fl  is  positive  definite  and  symmetric,  the  matrix  J/  -h  O  ' j  will 
also  possess  these  properties;  consequently,  it  has  a  square  root  decomposition  of  the  form 


{[/.Q-rf=[/.o-r. 


(75a) 


(75b) 


A  new  variable  of  integration  can  then  be  introduced  into  equation  (74)  via  the  formula 
z  =  1^7  +  Q“'  J'  ^  w.  Performing  the  necessary  algerbraic  manipulations  and  using  the  "Mixed 
Multiplication  Rule"  described  in  appendix  A  (see  equation  (A-8)}  yields 
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(76) 


A^^w;m,[/ +  n]  ^{u)du 


The  desired  integral  equation  for  y/  ^  (w)  follows  immediately  by  comparing  the  right-hand  side  of 
equation  (76)  with  equation  (69);  viz.. 


+  ^^y/ ^{v)du 


(77) 


=  |j^/  +  Q  ']  +  Q 


k  =  0,1,2 . 


The  result  embodied  in  equation  set  (77)  will  be  used  in  section  3  of  this  report  to  provide  a 
closed-form  solution  to  the  problem  of  detecting  random  multivariate  signals  embedded  in 
Gaussian  noise. 


2.4  DISCUSSION  OF  SECTION  2  RESULTS 

With  the  aid  of  the  /7-dimensional  Dirac  delta  function  representation  given  by  equation  (36), 
it  was  a  straightforward  task  to  obtain  the  Gram-Charlier  series  expansion  for  an  arbitrary 
multivariate  PDF  as  depicted  in  equation  (42).  Similarities  between  the  univariate  PDF  expansion 
and  its  multivariate  counterpart  are  clearly  evident.  In  the  univariate  case,  the  expansion 
coefficients  are  all  scalar  quantities;  the  first  three  values  being  1,0,  and  0,  respectively,  while  the 
remaining  coefficients  depend  successively  upon  progressively  higher-order  moments  starting  with 
third  order.^  An  examination  of  equation  set  (65)  reveals  that  the  same  pattern  prevails  in  the 
multivariate  case,  only  now  the  expansion  coefficients  are  vectors  of  exponentially  increasing 
dimension.  Furthermore,  using  only  the  first  term  of  the  Gram-Charlier  series  expansion  to 
approximate  an  actual  PDF,  whether  in  the  univariate  or  multivariate  case,  will  result  in  the 
equality  of  corresponding  moments  up  to  and  including  second  order.  As  might  be  expected, 
including  successively  more  terms  of  the  series  expansion  will  result  in  the  equality  of 
progressively  higher  order  moments. 


17/18 
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3.  MULTIVARIATE  RANDOM  SIGNAL  DETECTION  IN  GAUSSIAN  NOISE 


3.1  FORMULATION  OF  DETECTION  PROBLEM 

The  sonar  detection  problem  is  usually  formulated  as  a  binary  statistical  hypothesis  test.'^  ’* 
Under  the  null  hypothesis  Ho,  observations  (i.e.,  data  describing  the  output  of  acoustic  sensors) 
consist  of  noise  alone;  under  the  alternative  hypothesis  Hx,  observations  consist  of  noise  and 
signals  that  represent  responses  of  the  acoustic  sensors  to  the  presence  of  a  target  or  other 
interesting  object.  In  practice,  the  signals  and  noise  are  usually  time-varying  functions;  however, 
the  time  parameter  can  actually  be  suppressed  without  loss  of  generality  via  judicious  application 
of  the  Karhunen-Loeve  expansion.’^  Consequently,  only  the  single-time  epoch  problem  will  be 
considered,  the  results  being  readily  extended  to  time-varying  situations  without  undue 
complication.  With  this  in  mind,  the  aforementioned  detection  problem  is  mathematically 
described  as  follows: 


Ho-r^n, 

(78a) 

//,  :r  =  5+«, 

(78b) 

where 

(79) 

s  =  [s„s„...,s„]\ 

(80) 

«  =  k.«2. 

(81) 

Here  r,  denotes  the  output  of  the  sensor  element  (eg.,  the  response  of  the  /‘'’transducer  of  an 
acoustic  array),  and  5,  and  n,  are  the  corresponding  signal  and  noise  components,  respectively. 
The  vector  r  is  a  random  process  that  describes  the  output  (voltage)  of  an  array  of  acoustic 
transducers  comprising  the  data  The  signal  and  noise  vectors  are  assumed  to  be  uncorrelated 
random  processes,  with  the  noise  having  a  multivariate  Gaussian  PDF  given  by 

=  (82) 


The  signal  PDF,  denoted  by  Ps{^,  is  an  arbitrary  multivariate  distribution  that  satisfies  the 
relations 

1  (83a) 
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(83b) 


and,  additionally,  fulfills  all  the  usual  requirements^  of  a  legitimate  PDF. 

The  optimal  solution  to  the  binary  detection  problem  specified  by  equation  (78)  is  well- 
known  and  readily  available  in  the  literature.^’'’’'*'’'^  It  is  called  the  "likelihood  ratio"  test,  denoted 
by  A(r) ,  and  can  be  expressed  as 

When  both  signal  and  noise  vectors  are  Gaussian  random  processes,  the  integral  appearing  on  the 
right-hand  side  of  equation  (84)  can  be  evaluated  explicitly.  In  this  case,  a  simple  canonical 
representation  for  the  likelihood  ratio  is  obtained;  viz., 


A.fe) 


A't.d.r.) 


The  subscript  "g"  is  used  here  to  underscore  the  Gaussian  character  of  both  signal  and  noise.  The 
case  involving  non-Gaussian  signals  embedded  in  Gaussian  noise  yields  a  more  formidable 
expression  for  the  likelihood  ratio  that  is  not  amenable  to  closed-form  evaluation.  As  a  result, 
some  type  of  approximation  technique  is  required.  One  particularly  attractive  procedure  is 
described  in  the  next  section. 


3.2  DERIVATION  OF  A  CANONICAL  REPRESENTATION  FOR 
THE  LIKELIHOOD  RATIO 

The  derivation  of  a  canonical  representation  for  the  likelihood  ratio  A(r)  is  based  on  the 

result  embodied  in  equation  (84).  In  that  expression,  closed-form  evaluation  of  the  integral  is 
prohibited  because ps{s)  is  specified  as  an  arbitrary  PDF.  However,  this  difficulty  is  circumvented 
by  expanding ps{^  in  a  multivariate  Gram-Charlier  series  as  described  by  equation  (42);  viz., 

p/  s)  =  v(5,-  5,  r  (5,-  5,  r  (86) 

it=0  "  • 

where 
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and  are  the  (a?^  x  1  )  vector  valued  functions  defined  by  equation  set  (37).  If  equations 

(82)  and  (86)  are  substituted  into  equation  (84),  and  the  summation  and  integration  processes 
interchanged  (which  is  permissible  since  both  are  assumed  to  be  absolutely  convergent),  it  follows 
that 


k=0 


N{ 

(s:,s,L) 

m 

i 

1 

(88) 


Next,  observe  that  the  identity 


n( 

(sir,) 

N\ 

1 

-1/2 


(89) 


where 


(90) 


and 


(91) 


can  be  readily  deduced  via  straightforward  algebraic  manipulation  and  use  of  equation  (85).  If 
equations  (38)  and  (89)  are  now  substituted  into  the  integral  appearing  on  the  right-hand  side  of 

equation  (88),  and  the  variable  of  integration  changed  via  the  formula  m  =  '^^(^“1). 


A(/:)  =  A,(rE|7£ 

A:=0  ^  ” 


(92) 


The  integral  appearing  in  equation  (92)  is  identical  to  the  one  shown  in  equation  set  (77),  and 
the  matrix  Q  defined  by  equation  (91)  is  symmetric  and  positive  definite.  As  a  result,  equation 
(92)  may  be  rewritten  in  the  form 


^r) 


k=0  ^  ■ 


[/  +  Q  '] 


■1/2 


(93) 
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Alternatively, 


A::xO  f  -i 


(94) 


where  the  latter  expression  follows  by  substituting  equations  (90)  and  (9 1)  into  equation  (74), 
algebraically  manipulating  the  result,  and  then  using  equation  (70)  along  with  the  “Mixed 
Multiplication  Rule”  described  in  appendix  A.  This  is  the  desired  canonical  representation  for  the 
likelihood  ratio.  It  might  be  noted  that  the  first  two  terms  of  the  series  given  by  equation  (94)  are 
explicitly  evaluated  in  appendix  B  (see  equations  (B-20)  through  (B-25)). 


3.3  DISCUSSION  OF  SECTION  3  RESULTS 


It  is  noteworthy  that  the  infinite  series  appearing  in  equation  (94)  can  be  summed  explicitly 


in  the  absence  of  the  matrix  Kronecker  product 


1/2 


.  In  particular. 


substituting  equation  sets  (33)  and  (38)  into  equation  (86)  reveals  that 


fc=0  ■ 


r  +r. 


r“ 

r. 

A(l  +  r,"'[F,+Fj 

D-J 

F 

{l+«}.[r,+F„]j 

(95) 


This  result  suggests  some  interesting  possibilities  for  simply  approximating  A(r)  in  situations 

where  the  signal-to-noise  ratio  is  high.  For  example,  an  intuitively  pleasing  formula  is  obtained  by 
introducing  the  approximation  T  ^  +  T  „  «  T  ^ .  Under  this  assumption,  equations  (82),  (85), 

(94),  and  (95)  may  be  combined  and  subsequently  manipulated  to  produce  the  expression 


A(,). 


Psjr-^ 


(96) 


Comparing  equation  (84)  with  equation  (96)  reveals  that  the  noise  PDF  takes  on  Dirac  delta 
function  characteristics  as  the  signal  strength  increases  to  the  point  of  overwhelming  the  noise. 
This  is  expected,  since  the  condition  T  ^  +  T  „  F  ^  has  a  mathematically  analogous  effect  on  the 

integration  process  depicted  in  equation  (84)  as  the  condition  F  „  [O].  Recalling  equation  (82) 

and  the  fact  that 


rHo|  a({!;  -  s); ».  r  J  i  -  n), 


(97) 


readily  explains  the  approximation  given  by  equation  (96). 
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Equation  (95)  may  also  be  used  to  confirm  earlier  results  depicted  by  equation  (85). 
Indeed,  for  Gaussian  signal  and  noise  vectors,  the  right-hand  side  of  equation  (95)  reduces  to 
unity,  while  all  terms  on  the  left-hand  side  vanish  except  for  the  one  corresponding  to  A:  =  0.  In 
this  case,  the  equivalence  between  (85)  and  (94)  is  manifest.  Finally,  even  if  no  constraints  are 
imposed  on  the  signal-to-noise  ratio,  the  canonical  representation  depicted  in  equation  (94)  may 
be  recast  in  the  form 


A(r)  =  . 


r 

/’,(£  +  r/''[r,  +  r„p{r -£-«}) 


- (dZ  Irlr, ''’®-  •  I  &{.(!!.  {I + 2).[r,  +  r  J. 


where 


[r,+r„p®...' 

L  -  J 

(98) 


(99) 


In  practical  applications,  situations  may  arise  wherein  an  A^-term  series  approximation  based 
on  equation  (98)  will  provide  more  accuracy  than  a  similar  type  approximation  based  on  equation 
(94).  Of  course,  the  opposite  situation  may  also  arise  under  appropriate  circumstances. 
Ostensibly,  equation  (94)  should  prove  to  be  more  utilitarian  under  low  signal-to-noise  ratio 
conditions,  while  equation  (98)  is  preferred  when  high  signal-to-noise  ratio  conditions  prevail. 
Having  both  representations  available  thus  provides  a  greater  degree  of  flexibility  for  dealing  with 
specific  problems. 


4.  SUMMARY  AND  CONCLUSIONS 


Section  2  of  this  report  extended  the  well-known  Gram-Charlier  series  expansion  technique 
for  an  arbitrary  univariate  PDF  to  the  multivariate  case.  The  approach  employed  here  deliberately 
avoided  explicit  use  of  tensor  analysis  and  multivariate  Hermite  polynomials,  both  of  which  are 
cumbersome  to  apply,  even  under  the  best  of  circumstances.  The  desired  expansion  was  obtained 
by  developing  a  suitable  infinite  series  representation  for  the  ^-dimensional  Dirac  delta  function, 
and  then  judiciously  using  Kronecker  product  formulae  and  matrix  calculus  rules  to  manipulate 
the  result. 

The  analogy  between  a  univariate  PDF  series  expansion  and  its  multivariate  counterpart  is 
clearly  evident.  Expansion  coefficients  in  the  univariate  case  are  all  scalar  quantities,  the  first 
three  values  being  1,0,  and  0,  respectively.  The  remaining  coefficients  depend,  in  succession,  on 
progressively  higher  order  moments  starting  with  the  third  order.  For  the  multivariate  case,  a 
similar  pattern  exists;  however,  the  expansion  coefficients  are  now  vector  quantities  whose 
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dimension  increases  in  an  exponential  manner  with  the  summation  index.  In  both  cases,  only  a 
single  term  in  the  series  expansion  is  needed  to  exactly  replicate  any  Gaussian  PDF.  Under  such 
circumstances,  all  the  expansion  coefficients  vanish  except  the  first.  When  the  actual  PDF  is  non- 
Gaussian,  using  a  one-term  series  approximation  leads  to  the  equality  of  corresponding  moments 
up  to  and  including  second  order.  As  might  be  expected,  retaining  successively  more  terms  in  the 
series  approximation  of  the  actual  PDF  will  result  in  the  equality  of  progressively  higher  order 
moments. 

While  the  computational  procedure  for  explicitly  evaluating  the  multivariate  expansion 
coefficients  is  by  no  means  trivial,  it  is  nevertheless  tractable,  A  particularly  attractive  feature  of 
the  formulae  is  that  the  statistical  mean  and  covariance  associated  with  the  actual  PDF  always 
appear  in  their  compact  vector  and  matrix  formats,  respectively.  The  need  to  explicitly  depict 
individual  vector  components  or  matrix  elements  has  been  avoided  by  using  Kronecker  products. 
This  artifice  greatly  simplifies  expressions  that  would  otherwise  be  unwieldy  and  not  at  all 
amenable  to  analysis  and/or  interpretation. 

An  important  application  of  the  multivariate  Gram-Charlier  series  expansion  technique  is  in 
solving  the  binary  detection  problem  for  multivariate  signals  embedded  in  Gaussian  background 
noise.  Section  3  of  this  report  addressed  that  problem  in  the  context  of  sonar  array  signal 
processing.  The  solution  was  specified  as  a  binary  statistical  hypothesis  test,  viz.,  the  "likelihood 
ratio  test."  A  canonical  representation  for  the  likelihood  ratio  was  developed  in  the  form  of  an 
infinite  series  whose  individual  terms  depend  on  the  received  measurement  vector  together  with 
the  signal  and  noise  statistics.  A  closed-form  expression  for  each  term  was  subsequently  obtained 
via  use  of  an  integral  equation  satisfied  by  the  multivariate  Gram-Charlier  expression  coefficients. 
This  description  reveals  manifest  similarities  between  the  multivariate  and  univariate  canonical 
representations.  In  fact,  it  is  easily  demonstrated  that  the  former  reduces  to  the  latter  as  the 
dimensionality  parameter  n  approaches  unity. 

Another  important  consequence  of  evaluating  the  individual  terms  of  the  infinite  series  in 
closed-form  is  that  it  facilitates  the  derivation  of  likelihood  ratio  approximations  that  are  suitable 
for  use  under  various  operating  conditions  In  particular,  the  aforementioned  description  of 
individual  terms  allows  the  original  infinite  series  to  be  decomposed  into  two  other  infinite  series, 
one  of  which  is  summable  to  a  recognizable  function.  As  a  result,  the  likelihood  ratio  can  be 
represented  in  either  of  two  canonical  forms.  The  forms  are  equivalent  and  of  identical 
composition;  viz.,  both  are  expressed  as  infinite  series.  Ostensibly,  the  convergence 
characteristics  of  these  two  series  are  quite  different.  When  the  signal-to-noise  ratio  is  high,  one 
series  appears  to  converge  rapidly,  while  the  other  does  not.  When  the  signal-to-noise  ratio  is 
low,  this  pattern  seems  to  reverse  itself  Consequently,  in  practical  applications,  judiciously 
choosing  the  appropriate  canonical  representation  before  making  an  A^term  series  approximation 
should  yield  a  more  accurate  description  of  the  actual  likelihood  ratio.  Having  both 
representations  available  obviously  enhances  the  capability  to  address  specific  problems  more 
effectively. 
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APPENDIX  A 


KRONECKER  PRODUCT  FORMULAE  AND  MATRIX  CALCULUS  RULES 


Presented  in  this  appendix  are  some  formulae  for  manipulating  Kronecker  products  and  for 
performing  matrix  calculus.  The  results  are  presented  without  proof  Interested  readers  are 
referred  to  Graham  for  detailed  derivations  and  a  comprehensive  treatment  of  the  subject 
material. 


DEFINITION  OF  THE  KRONECKER  PRODUCT 


Consider  a  matrix  A  =  of  order  (m  x  p)  and  a  matrix  B  =  [bij\  of  order  (r  x  s).  The 
Kronecker  product  of  the  two  matrices,  denoted  by  A  ®  B,  is  defined  as  the  partitioned  matrix: 


a,, 5 

a,,B  . 

a^.B  . 

..  a^^B 

A®B^ 


a^,B  a_^B 


m2'‘ 


mp 


(A-1) 


A  0  B  is  seen  to  be  a  matrix  of  order  imr  x  ps).  It  has  mp  blocks,  the  block  is  the  matrix 
a,j  B  of  order  (rxs). 

Some  properties  and  rules  for  Kronecker  products  are: 


1 .  (aA)<8>  (J3B)  =  a^A  0  B)  for  any  scalars  a,  p.  (A-2) 

2.  If  A  and  B  are  both  (m  x  p),  and  C  is  (r  x  s),  then 

iA+B)®C  =  A®C  +  B®C.  (A-3) 

3.  If  A  is  (mxp)  and  both  B  and  C  are  (r  x  5),  then 

A<S)(B  +  C)  =  A®B+A®C.  (A-4) 

4.  A®(B0Q  =  (A®B)0C.  (A-5) 

5.  (A0Bf  =  A^0B'^.  (A-6) 

6.  If  A  and  B  are  both  (m  x  m),  then 

tr(A  0B)  =  tr{A)  tr{B),  (A-7) 
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where  tr{  )  is  the  trace  operator. 

7.  If  is  (tw  X p\  B  \s.{qx  r),  C  \&{px  5),  and  Z)  is  (r  x  v),  then 
{A  ®  B){C  ®D)  =  (AQ®  (BD). 

8.  A  and  B  are  both  nonsingular,  then 
(A®  By'  =  A-'®B-'. 

9.  If  A  is  (mx  p)  and  .5  is  x  r),  then 
A®B  =  U(m,q)  (B  ®  A)  U{r,py 

where 

=  i  Iff, (*)/(/)}»{£, ('KW) 

/-I  >=] 

is  the  (Jcl  X  kJ)  permutation  matrix*^  defined  for  all  integer  values  of  k  and  /,  and 


are  («  x  1)  unit  vectors  defined  for  all  integer  values  of  n 
Let 


d  _  d  d  d 
du  du^  ’  du^  ’  ’  du^ 


for  u-{^,x), 


(A-8) 
(A-9) 
(A- 10) 

(A-11) 

(A-12) 


(A-13) 


and  assume  that  F  is  a  symmetric,  positive  definite  {n  x  n)  matrix.  Under  these  assumptions,  it  is 
easily  shown'®  that  F  has  the  square  root  decomposition 

F=:F'^'F''^  (A-14) 


A-2 


(A-15) 


Now  observe  that 


d 

—  I  ri 


\dx) 


(  S 

pl/2  ^ 

K  ffxJ 


Vr 


-1/2 


dx) 


(A- 16) 


Consequently,  the  square  of  this  scalar  differential  operator  can  be  expressed  in  the  form 


\dx) 
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n2 


Kdx) 


(  ^Vr 


'1/2 
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1/2 

V  dx) 


0 


A 


-1/2 


dx) 


.1/2  d 


dx) 


(A-17) 


since  the  Kronecker  product  of  two  scalar  quantities  is  identical  to  the  ordinary  scalar  product. 
Using  equation  (A-8)  allows  equation  (A-17)  to  be  rewritten  as 


J 


dx) 


Kdx) 


dx) 


d^ 

pl/2  ^ 

<  dx) 


-|2 


(A-18) 


where  the  notation 


[m..f  =  A®  A®...®  A 


^  =  0,1,2,. 


k  times 


has  been  employed.  Repeating  this  process  yields  the  following  result  by  induction: 


(A- 19) 


^  d^ 


Vdx) 
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\dx) 


dV 

d^J 


k  r 


1/2  ^ 
dxj 


0... 


(A-20) 


vec{  )  OPERATOR 
If 

^=[ai  :_£2  :  (A-21) 

is  an  (/n  x  p)  matrix  partitioned  as  shown,  then 
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«1 


G.2 


vec[A)  = 


(A-22) 


is  a  {mp  X  1)  vector.  Some  useful  formulae  connecting  the  vec  of  a  matrix  product  with  the 
Kronecker  product  and  the  trace  of  a  matrix  product  are  specified  below: 

1 .  If  5  is  (/w  X  n),  C  is  («  X  q),  and  D  is  x  5),  then 


vec{BCD)  =  (D^  ®  B)  vec{C).  (A-23) 

2.  If  5  is  (w  X  «)  and  m  is  a  («  x  1)  vector,  then 

vec{B  <S>u}  =  vec{B)  0  u  (A-24) 

vec{ui/)  =  u'S>  u.  (A-25) 

3 ,  If  5  is  (w  X  rt)  and  C  is  («  x  m),  then 

{vec(B^)y  {vec(C)}  =  tr(BC).  (A-26) 


DIFFERENTIATION  FORMULAE  FOR  KRONECKER  PRODUCTS 

1 .  If  A  is  (m  xp).  Bis  (p\  q),  and  m  is  a  («  x  1 )  vector,  then 


B  +  {I®A][—®B 


(A-27) 


2.  If  A  is  {m  X  p),  B\s{qx  r),  and  w  is  a  («  x  1 )  vector,  then 


-y®{A®B] 


,  -  0  ^ 

\du 


05  +  |/  0  f/(/w,<7)) 


(A-28) 


In  equations  (A-27)  and  (A-28),  1  represents  the  (n  x  n)  identity  matrix,  while  U{m,q)  and  U{r,p) 
represent  permutation  matrices  defined  by  equation  (A-1 1),  Some  special  results  are  given  as 
follows: 
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If  >4  is  a  constant  (m  x  ri)  matrix  and  m  is  a  («  x  1)  vector,  then 


vec{u)  =  u , 

(A-29) 

vec(u^  )  =  u, 

(A-30) 

vec(t/A^  )-Au, 

(A-31) 

d  T 

—  ®u=I, 
du 

(A-32) 

—  ®u  =  vedlX 
du  -  ^  ’ 

(A-33) 

—  0  =  2m, 

du  - 

(A-34) 

—  0  {Au^  =  vec^A), 

(A-35) 

—  ®Ua^]  =  a\ 

du  > 

(A-36) 

If  ^  is  a  constant  («  x  n)  matrix  and  w  is  a  («  x  1)  vector,  then 

-^®^uJA^  =  Au  +  A^u.  (A-37) 
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APPENDIX  B 

COMPUTATION  OF  THE  FIRST  TWO  NON-TRIVIAL  TERMS  OF  THE 
MULTIVARIATE  GRAM-CHARLIER  SERIES 


In  this  appendix,  the  first  two  non-trivial  terms  appearing  in  the  Gram-Charlier  series  given 
by  equation  (42)  are  evaluated.  Since  each  term  in  that  series  takes  the  form  of  a  scalar  product 
of  two  vectors,  various  mathematical  artifices  may  be  employed  to  simplify  the  calculations.  To 
this  end,  observe  that  the  first  two  non-trivial  terms  are  associated  with  the  indices  A:  =  3  and  k  = 
4.  Those  corresponding  to  A:  =  0,  A  =1,  and  k  =  2  art  trivial  to  compute;  their  numerical  values 
being  1,  0,  and  0,  respectively.  Focusing  on  the  non-trivial  term  with  index  A  =  3,  the  evaluation 
process  begins  by  using  the  formula 

y/ ^{u)  =  u'S) u~ 

which  follows  from  equations  (57d),  (58c),  and  (A-25).  Substituting  equation  (B-1)  into 
equation  (58d),  and  combining  that  result  with  equations  (57d),  (A-23),  and  (A-24)  yields 

y/^{u)  =  {/0C}|vec(/)®  -I-M0  vec(/)- M®  w®  M.  (B-2) 

Aided  by  equations  (B-2)  and  (A-8),  along  with  the  relations 

^3{/®C}{vec(/)®  w)  =l\{2vec{l)®u],  (B-3) 


and 


/?^{m®  vec(/)|  =y5^{vec(/)®  m|,  (B-4) 

it  can  be  shown  that 

=  |£  [(x -x)®...j^p(x)t/x|  ||z® z- 3vec(r  '^J® z|,  (B-5) 

where 

z  =  r"’{x-x}.  (B-6) 

From  equations  (A-24)  and  (A-25)  it  also  follows  that 

|r®z-3vec(r"'jj®z  =  vecjjzz^  - 3r  ' j®z|,  (B-7) 


and 
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[(x  -  x)®. . .]  =  vec(Ay, 


(B-8) 


where 


^3 


(B-9) 


is  the  matrix  of  third  order  central  moments  associated  with  /7(x).  Finally,  substituting  equations 
(B-7)  and  (B-8)  into  equation  (B-5),  and  using  equation  set  (38)  and  equation  (A-26)  leads  to  the 
desired  result;  viz., 


(B-IO) 


Computation  of  the  non-trivial  term  associated  with  index  /:  =  4  is  initiated  by  using  the 
formula 


(B-11) 


in  conjunction  with  equation  (58e)  to  obtain 


(B-12) 


Combining  this  expression  with  equation  set  (63),  and  then  using  equations  (A-8)  and  (A-23) 
along  with  some  algebraic  manipulations,  yields  the  result 


^,(x;x,r)  = 


—  ®— 
2  2 


[r'"'®...]%ec{Q(z)). 


(B-13) 


where 


Q(z)  =  |v^^^  -r-')|{vec(^''  -r-')}^  -4r-’  ®(^^)  +  2r-’  ®r-’.  (b-m) 

and  z  is  defined  by  equation  (B-6).  With  the  aid  of  equations  (B-13)  and  (65e),  and  the  identity 


4f  ®f 


(B-15) 


it  follows  that 
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(B-16) 


f,  {x;  r)  =  [|^^  [{x  -  x)0. . .]'  p{x)dx  -  vec{  c(r  0  r)) 

-  vec(r)0vec(r)]  vec  {«<?))■ 


Next,  observe  the  relation 


|vec  {c(r®r)}f  vec  {fXd)  =  2[vec(r )  0  vec{r  )  vec{!D.(z)]. 


(B-17) 


A4  =  [{x  -  x)(x  -  f )" ]  0  [(x  -  x)(x  -  xf  ]/<x  y X 


(B-18) 


is  defined  as  the  matrix  of  fourth  order  central  moments  associated  with  p(x),  then  equations 
(B-17),  (B-18),  and  (A-26)  allow  equation  (B-16)  to  be  rewritten  in  the  desired  form;  i.e., 


=  /r|[A,  -3{vec(r)}{vec(r)}"]Q(z)|. 


(B-19) 


Equations  (B-10)  and  (B-19)  represent  simplified  expressions  for  the  first  two  non-trivial 
terms  in  the  Gram-Charlier  series  expansion  of p(x).  It  is  noteworthy  that  analogous  expressions 
also  exist  for  the  corresponding  non-trivial  terms  in  the  likelihood  ratio  representation  given  by 
equation  (94).  They  are 


2T(r"[r,+r.p)®...J|,(r,||+2),[r,+r,]) 


(r|A-,rU"-3[r,+r,]'')»yl|, 


(B-20) 


^ rh ■  •]  ^4 [r. + r, ]) 

=  tr\  A;-3{vec(r^)}{vcc(r,)}nn*(^) 


(B-21) 


where 


(B-22) 


B-3 


(B-23) 


a;=  £ j(5  - 1)(5  - 1)"]  ®  (^  -  if  PsM)^^’ 

^4=^  (B-24) 

n*(^)-  |vec|>y^  -[r,  +  rj  -[r,  +  rj  '|| 

(B-25) 

-  4[r,  +  ]  ’  ®  (^") + 2[r^ + r J'  ®  [r^  +  r J’ . 
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